#!/usr/bin/env julia
"""
NSEMSolver.jl Test Suite

Comprehensive test suite for the NSEMSolver.jl package covering all major
functionality including manufactured solutions, adaptive refinement, and
solver backends.
"""

using Test
using LinearAlgebra
using NSEMSolver

# Check if we should run only Aqua tests
run_aqua_only = length(ARGS) > 0 && ARGS[1] == "aqua"

@testset "NSEMSolver.jl Tests" begin
    
    if run_aqua_only
        # Run only Aqua quality tests
        include("test_aqua.jl")
        return
    end
    
    @testset "Package Loading and Initialization" begin
        @test isdefined(NSEMSolver, :solve_navier_stokes_2d)
        @test isdefined(NSEMSolver, :NSOptions)
        @test isdefined(NSEMSolver, :NSResult)
        @test isdefined(NSEMSolver, :create_manufactured_solution)
        
        # Test thermal-fluid physics exports
        @test isdefined(NSEMSolver, :BoussinesqOptions)
        @test isdefined(NSEMSolver, :HeatTransferOptions)
        @test isdefined(NSEMSolver, :MultiPhysicsResult)
        @test isdefined(NSEMSolver, :solve_coupled_flow_heat_2d)
        @test isdefined(NSEMSolver, :ThermalBoundarySystem)
        @test isdefined(NSEMSolver, :SmagorinskyModel)
        
        # Test backend detection
        @test isdefined(NSEMSolver, :HAS_PETSC)
        @test isdefined(NSEMSolver, :HAS_GCR)
        @test isdefined(NSEMSolver, :HAS_MPI)
        
        println("✅ Package components loaded successfully")
    end
    
    @testset "NSOptions Configuration" begin
        # Test default options
        options_default = NSOptions()
        @test options_default.N >= 2
        @test options_default.nu > 0
        @test options_default.tfinal > 0
        @test options_default.cfl > 0
        
        # Test custom options
        options_custom = NSOptions(
            N = 5,
            nu = 0.02,
            tfinal = 2.0,
            solver = :julia,
            verbose = false
        )
        @test options_custom.N == 5
        @test options_custom.nu == 0.02
        @test options_custom.tfinal == 2.0
        @test options_custom.solver == :julia
        @test options_custom.verbose == false
        
        println("✅ NSOptions configuration tests passed")
    end
    
    @testset "Manufactured Solutions" begin
        # Test 2D manufactured solution
        @testset "2D Manufactured Solution" begin
            sol_2d = create_manufactured_solution(2)
            @test sol_2d.dim == 2
            @test isa(sol_2d, NSManufacturedSolution)
            
            # Test evaluation at a point
            x, y, z, t = 0.5, 0.3, 0.0, 0.5
            u_val = calculate_manufactured_u(sol_2d, x, y, z, t)
            v_val = calculate_manufactured_v(sol_2d, x, y, z, t)
            w_val = calculate_manufactured_w(sol_2d, x, y, z, t)
            p_val = calculate_manufactured_p(sol_2d, x, y, z, t)
            
            @test isa(u_val, Real)
            @test isa(v_val, Real)
            @test w_val == 0.0  # Should be zero for 2D
            @test isa(p_val, Real)
            
            # Test forcing terms
            fx, fy, fz = calculate_manufactured_forcing(sol_2d, x, y, z, t, 0.01)
            @test isa(fx, Real)
            @test isa(fy, Real)
            @test fz == 0.0  # Should be zero for 2D
            
            println("✅ 2D manufactured solution tests passed")
        end
        
        # Test 3D manufactured solution
        @testset "3D Manufactured Solution" begin
            sol_3d = create_manufactured_solution(3)
            @test sol_3d.dim == 3
            
            # Test evaluation at a point
            x, y, z, t = 0.4, 0.6, 0.2, 0.3
            u_val = calculate_manufactured_u(sol_3d, x, y, z, t)
            v_val = calculate_manufactured_v(sol_3d, x, y, z, t)
            w_val = calculate_manufactured_w(sol_3d, x, y, z, t)
            p_val = calculate_manufactured_p(sol_3d, x, y, z, t)
            
            @test isa(u_val, Real)
            @test isa(v_val, Real)
            @test isa(w_val, Real)  # Should be non-zero for 3D
            @test isa(p_val, Real)
            
            println("✅ 3D manufactured solution tests passed")
        end
        
        # Test divergence-free condition
        @testset "Divergence-Free Validation" begin
            sol_2d = create_manufactured_solution(2)
            is_div_free_2d = validate_manufactured_solution(sol_2d)
            @test is_div_free_2d
            
            sol_3d = create_manufactured_solution(3)
            is_div_free_3d = validate_manufactured_solution(sol_3d)
            @test is_div_free_3d
            
            println("✅ Divergence-free validation tests passed")
        end
    end
    
    @testset "SEM Operators" begin
        # Test SEM operator creation
        sem_ops = create_sem_operators(4, 2, 1.0, 2)
        
        @test isa(sem_ops, SEMOperators)
        @test length(sem_ops.nodes) > 0
        @test size(sem_ops.D1, 1) == size(sem_ops.D1, 2)
        @test size(sem_ops.D2, 1) == size(sem_ops.D2, 2)
        @test length(sem_ops.coordinates) == 2  # 2D
        @test length(sem_ops.grids) == 2  # 2D
        
        println("✅ SEM operators tests passed")
    end
    
    @testset "Multi-Domain Geometry" begin
        # Test multi-domain creation
        multidomain = create_multidomain(3, 4, 0.01, 2; adaptive_refinement=true)
        
        @test isa(multidomain, MultiDomain)
        @test multidomain.dim == 2
        @test multidomain.n == 4
        @test multidomain.index_sum > 0
        @test size(multidomain.matrix_dist) == (5, 5)  # n_block + 2
        
        # Test refinement capabilities
        @test multidomain.adaptive_refinement == true
        @test multidomain.refinement_map !== nothing
        @test length(multidomain.corner_points) > 0
        
        # Test without adaptive refinement
        multidomain_no_ref = create_multidomain(2, 3, 0.01, 2; adaptive_refinement=false)
        @test multidomain_no_ref.adaptive_refinement == false
        @test multidomain_no_ref.refinement_map === nothing
        
        println("✅ Multi-domain geometry tests passed")
    end
    
    @testset "L-Shaped Domain and Refinement" begin
        # Test L-shaped domain definition
        matrix_dist, corner_points = define_lshape_domain(3, 2)
        
        @test size(matrix_dist) == (5, 5)  # n_block + 2
        @test length(corner_points) > 0
        @test all(length(cp) == 2 for cp in corner_points)  # 2D corners
        
        # Test refinement map creation
        refinement_map = create_refinement_map(matrix_dist, corner_points, 3, 2)
        @test size(refinement_map) == size(matrix_dist)
        @test maximum(refinement_map) <= 3
        
        # Test refinement smoothing
        smooth_map = smooth_refinement_transitions!(copy(refinement_map))
        @test size(smooth_map) == size(refinement_map)
        
        println("✅ L-shaped domain and refinement tests passed")
    end
    
    @testset "Linear Solvers" begin
        # Create a small test matrix (use identity matrix for perfect conditioning)
        n = 10
        A = Matrix{Float64}(I, n, n)  # Identity matrix for perfect conditioning
        b = rand(n)
        x0 = zeros(n)
        
        # Test Julia solver creation
        julia_solver = create_ns_solver(A, :pressure, :julia, false)
        @test isa(julia_solver, NSLinearSolver)
        
        # Test solving with Julia backend
        x_sol = solve_linear_system!(julia_solver, b, x0, :pressure, true)  # Enable verbose
        @test size(x_sol) == size(b)
        residual_norm = norm(A * x_sol - b)
        println("Linear solver test - Residual norm: ", residual_norm)
        println("Linear solver test - Expected solution: ", b)
        println("Linear solver test - Actual solution: ", x_sol)
        @test residual_norm < 1e-6  # Check solution accuracy (relaxed for iterative solvers)
        
        # Test PETSc solver if available
        if NSEMSolver.HAS_PETSC
            petsc_solver = create_ns_solver(A, :pressure, :petsc, false)
            @test isa(petsc_solver, NSLinearSolver)
            
            x_petsc = solve_linear_system!(petsc_solver, b, x0, :pressure, false)
            @test size(x_petsc) == size(b)
            @test norm(A * x_petsc - b) < 1e-10
            
            println("✅ PETSc solver tests passed")
        end
        
        println("✅ Linear solver tests passed")
    end
    
    @testset "Small-Scale 2D Solver" begin
        # Test with minimal configuration for speed
        options = NSOptions(
            N = 3,
            n_block = 2,
            nu = 0.1,   # Higher viscosity for stability
            tfinal = 0.1,  # Short time
            cfl = 0.3,
            tol = 1e-4,  # Relaxed tolerance
            solver = :julia,
            adaptive_refinement = false,  # Disable for speed
            verbose = false
        )
        
        # Run solver
        result = solve_navier_stokes_2d(options)
        
        # Basic result validation
        @test isa(result, NSResult)
        @test result.u !== nothing
        @test result.v !== nothing
        @test result.p !== nothing
        @test result.solve_time > 0
        @test result.iterations >= 0
        
        # Check solution properties
        if result.converged
            @test maximum(abs.(result.u)) < 10.0  # Reasonable velocity magnitude
            @test maximum(abs.(result.v)) < 10.0
            @test maximum(abs.(result.p)) < 100.0  # Reasonable pressure magnitude
            println("✅ Small-scale 2D solver converged successfully")
        else
            @test result.residual_norm > 0
            println("⚠️  Small-scale 2D solver did not converge (expected for short time)")
        end
        
        println("✅ 2D solver basic functionality tests passed")
    end
    
    @testset "Error Indicators" begin
        # Create test velocity fields
        n = 5
        u_test = rand(n, n)  # Random u-velocity field (n x n)
        v_test = rand(n, n)  # Random v-velocity field (n x n)
        D1_test = rand(n, n) # Random derivative matrix
        D2_test = rand(n, n)
        
        # Test 2D error indicators
        error_2d = calculate_physics_error_indicators(u_test, v_test, nothing, D1_test, D2_test, 0.01, 2)
        @test size(error_2d) == size(u_test)
        @test all(error_2d .>= 0)  # Error indicators should be non-negative
        
        # Test 3D error indicators
        w_test = rand(n, n)  # Random w-velocity field (n x n)
        error_3d = calculate_physics_error_indicators(u_test, v_test, w_test, D1_test, D2_test, 0.01, 3)
        @test size(error_3d) == size(u_test)
        @test all(error_3d .>= 0)
        
        println("✅ Error indicator tests passed")
    end
    
    @testset "Utility Functions" begin
        # Test polynomial order mapping
        p_orders = [get_polynomial_order_from_refinement(i) for i in 0:5]
        @test all(p_orders .>= 2)  # Minimum polynomial order
        @test all(diff(p_orders) .>= 0)  # Should be non-decreasing
        
        # Test stable timestep calculation  
        sem_ops = create_sem_operators(3, 2, 1.0, 2)
        dt = compute_stable_timestep(sem_ops, 0.5, 0.01)
        @test dt > 0
        @test dt < 1.0  # Should be reasonably small
        
        println("✅ Utility function tests passed")
    end
    
    # Include comprehensive visualization tests
    println("\n🎨 Running Visualization System Tests...")
    include("test_visualization.jl")
    
    # Include thermal-fluid physics tests  
    println("\n🌡️  Running Thermal-Fluid Physics Tests...")
    include("test_thermal_physics.jl")
    
    # Include advanced solver tests
    println("\n🚀 Running Advanced Solver Tests...")
    include("test_advanced_solvers.jl")
    
    # Include package quality tests
    println("\n🔍 Running Package Quality Tests...")
    include("test_aqua.jl")
end

println("\n🎉 All NSEMSolver.jl tests completed!")
println("✅ Package is ready for use")